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Abstract 

We  report  serious  bugs  in  floating-point  computations  for  evaluating  elementary  functions  in  the 
Embedded  GNU  C  Library.  For  instance,  the  sine  function  can  return  values  larger  than  1053  in 
certain  rounding  modes.  Further  investigation  also  exposed  faulty  implementations  in  the  most  re¬ 
cent  version  of  the  library,  which  seemingly  fixed  some  bugs,  but  only  by  discarding  user-specified 
rounding-mode  requirements. 
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1  Introduction 


We  have  found  floating-point  bugs  in  Linux  systems  using  Embedded  GLIBC  (EGLIBC)  version 
2.16  or  older.  EGLIBC  is  a  variant  of  the  GNU  C  Library  (GLIBC)  which  is  used  as  the  default 
implementation  in  many  distributions  including  Debian,  Ubuntu,  and  their  variants. 

The  following  C  program  computes  the  value  of  sin(— 2.437592)  in  double-precision  after 
setting  the  rounding  direction  to  upward  (+oo). 


1 

♦include  <math.h> 

2 

♦include  <fenv.h> 

3 

♦include  <stdio.h> 

4 

5 

int  main ( )  { 

6 

double  x  =  -2.437592; 

7 

fesetround (FE_UPWARD) ; 

8 

print f (" sin (%f) =%f\n" ,  x,  sin (x) ) ; 

9 

return  0; 

10 

} 

The  IEEE754  standard  [1]  does  not  specify  correct  rounding  methods  on  elementary  func¬ 
tions  such  as  the  exponential,  logarithm,  and  trigonometric  functions.  Programmers  and  engi¬ 
neers  usually  expect  the  program  to  print  out  an  approximated  value  around  sin(— 2.437592)  ~ 
—0.64727239229  with  an  “acceptable”  amount  error.  However,  they  all  should  agree  that  the  result 
be  in  the  range  between  -1  and  1,  even  in  the  worst  case. 

However,  a  surprising  result  appears  if  we  compile  and  execute  the  program  in  a  machine 
running  Ubuntu  12.04  LTS  (or  any  system  with  EGLIBC-2.15).  The  value  is  greater  than  1053  and 
it  should  not  be  a  return  value  of  sine  function  in  any  sense. 


$  gcc  exp_bug.c  -lm  &&  ./a. out 

sin (-2. 4375 92) =191561981424 93 6 94 3059347 92 7032 14803 02 8731397 92 094 16704. 00000 


Here  is  another  C  program  computing  cosh(3. 113408)  with  directed  rounding  toward  Too. 
This  example  is  more  interesting  because  it  shows  different  results  on  Intel  and  AMD  machines, 
and  both  of  the  results  have  serious  problems. 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 


#include  <math.h> 
#include  <fenv.h> 
#include  <stdio.h> 


int  main ( )  { 

double  x  =  3.113408; 
f eset round (FE_UPWARD)  ; 

print f ( "cosh (%f)  =  %f\n",  x,  cosh (x) ) ; 
return  0; 

} 
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Table  1:  Experiment  Setup 


Function 

Domain 

Range 

Function 

Domain 

Range 

sin 

'  1  ' 

o 

co 

o 

C5 

o 

co 

o 

05 

[-1.0, 1.0] 

acos 

[-1.0, 1.0] 

[  oo,  +cx>] 

cos 

[_103°6,  10306] 

[-1.0, 1.0] 

asin 

[—1.0, 1.0] 

[  oo,  +cx>] 

tan 

[-10306,  10306] 

[-ex),  +oo] 

atan 

[-1.0, 1.0] 

[  oo,  +cx>] 

cosh 

[-500,  500] 

[1.0,  +oo] 

exp 

[-100,100] 

[0.0,  +CX)] 

sinh 

[-500,500] 

[-(X),  +cx>] 

log 

[l0-3°6,  10306] 

[  oo,  +cx>] 

tanh 

sqrt 

[-100,100] 
[0.0, 10306] 

[-1.0, 1.0] 
[0.0,  +oo] 

log  10 

[10~306, 10306] 

[  oo,  +cx>] 

In  a  machine  with  Intel  Core  i7  CPU,  the  program  outputs  inf  while  a  machine  with  AMD 
Opteron  processor  produces  —160.191709. 


[INTEL]  $  gcc  cosh_bug.c  -lm  &&  ./a. out 
cosh (3 . 113408)  =  inf 

[AMD]  $  gcc  cosh_bug.c  -lm  &&  ./a. out 
cosh(3. 113408)  =  -160.191709 


Note  that  cosh(3. 113408)  ~  11.2710174432  and  the  both  results  inf  and  —160.191709  are 
simply  wrong.  Moreover,  each  of  the  wrong  results  have  significant  implications: 

•  Intel  (inf):  It  has  a  contagious  effect  in  subsequent  computations,  inf  is  a  special  value  in  the 
IEEE754  standard  which  indicates  an  overflow  in  a  computation.  If  one  of  subexpressions  is 
evaluated  to  inf  then  in  general  the  main  expression  also  becomes  infinity  (+oo  or  —  oo)  or 
NaN  (Not  a  Number). 

•  AMD  (-160.191709):  Mathematically,  cosh(x)  is  greater  than  or  equal  to  1  for  all  x  €  M. 
As  a  result,  programmers  and  engineers  write  algorithms  based  on  the  invariant  \/x.  1  < 
cosh(x).  This  result,  —160.191709,  breaks  the  invariant  and  could  cause  an  unexpected 
behavior. 


2  Floating-point  Bugs  in  EGLIBC  (<  2.16) 

We  have  tested  the  following  math  functions  in  C  standard  library: 

sin,  cos,  tan,  acos,  as  in,  atan,  cosh,  sinh,  tanh,  exp,  log,  loglO,  sqrt. 

For  each  function  /,  we  take  100, 000  random  numbers  from  a  subset  of  function  /’ s  domain. 
Table  1  shows  each  function’s  sampling  domain  and  range.  We  pick  the  sampling  domain  carefully 
so  that  the  result  of  the  computation  can  be  fit  in  a  double-precision  variable.  We  consider  the  four 
rounding  modes  supported  by  C99  standard  [2]: 

•  (nearest),  -*0  (toward  zero),  f  (toward  +oo),  f  (toward  —  oo). 
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Table  2:  Experimental  results  on  Intel  and  AMD  machines:  f^,  f3,  and  f  *°  indicate  a  function  f 
with  rounding  mode  toward  +oo,  toward  —  oo,  and  toward  0  respectively.  “Inconsistent”  denotes 
the  number  of  cases  in  which  the  difference  of  two  results  are  larger  than  220  ULP  (Unit  of  Least 
Precision).  “Incorrect”  denotes  the  number  of  cases  in  which  ic(x)  is  out  of  /’ s  mathematical 


range.  “±oc 

Function 

)”  denotes  the  number  of  cases  in  which  fc 
Intel 

(x)  is  either  —  oo  or  +oo. 

|  AMD 

Inconsistent 

Incorrect 

±oo 

Total 

Inconsistent 

Incorrect 

±oo 

Total 

•  -t 

sin1 

10055 

446 

0 

10501 

9853 

450 

0 

10303 

sin'1' 

9619 

497 

0 

10116 

10009 

450 

0 

10459 

sin-*'0 

10087 

436 

0 

10523 

9904 

423 

0 

10327 

cos1' 

10097 

434 

0 

10531 

9880 

423 

0 

10303 

COS1' 

9815 

442 

0 

10257 

9910 

461 

0 

10371 

cos-*0 

9737 

444 

0 

10181 

9913 

441 

0 

10354 

tan1 

12218 

0 

0 

12218 

12452 

0 

0 

12452 

tan1' 

12387 

0 

0 

12387 

12378 

0 

0 

12378 

O 

T 

£ 

12486 

0 

0 

12486 

12506 

0 

0 

12506 

cosh1^ 

18768 

30139 

935 

49842 

37091 

12295 

291 

49677 

cosh1. 

49766 

0 

0 

49766 

49673 

0 

0 

49673 

cosh^0 

49713 

0 

0 

49713 

49772 

0 

0 

49772 

sirih1 

47807 

0 

0 

47807 

47451 

0 

266 

47717 

sinh-i 

47493 

0 

0 

47493 

47676 

0 

0 

47676 

sinh-*0 

47911 

0 

0 

47911 

48046 

0 

0 

48046 

tanh1 

3107 

0 

0 

3107 

3242 

0 

0 

3242 

tarih1 

3135 

0 

0 

3135 

3268 

0 

0 

3268 

exp1" 

47386 

2536 

0 

49922 

37883 

11708 

124 

49715 

exp1. 

49646 

0 

0 

49646 

49978 

0 

0 

49978 

exp^° 

50022 

0 

0 

50022 

49836 

0 

0 

49836 

For  each  sample  x  and  for  each  rounding  mode  rnd,  we  compute  two  values  /™d(a;)  and 
.Atpfr.  (a‘)  where  fc  is  a  function  /  in  C  standard  library  and  /mpfr  is  a  function  /  in  the  GNU 
MPFR  library.  MPFR  supports  arbitrary-precision  floating-point  computation  and  we  use  it  as  a 
reference  implementation  to  have  a  comparison.  The  correctness  of  MPFR  is,  of  course,  another 
issue  and  we  do  not  discuss  it  here.  In  the  experiments,  we  use  256-bit  precision  for  MPFR. 

We  have  the  following  expectations  for  the  two  values: 

•  Consistency:  The  difference  between  f™d(x)  and  /mpfr(x)  should  not  be  too  large.  In  this 
experiment,  we  set  the  threshold  of  220  ULP  (Unit  of  Least  Precision)  which  is  the  spacing 
between  two  adjacent  floating-point  numbers.  Note  that  IEEE754  double-precision  format 
has  53  bits  of  precision  and  220  ULP  implies  that  it  loses  20-bit  precision  out  of  53.  If 
|/™d(x)  —  /mpfr^)!  >  220ULP,  we  label  the  case  as  “inconsistent”. 

•  Correctness:  The  value  of  /™d  (x)  should  be  in  the  range  of  the  mathematical  function  /. 
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For  instance,  sin™d(a;)  has  to  be  between  -1.0  and  1.0  no  matter  how  imprecise  it  is. 

We  run  the  experiments1  on  two  machines  -  one  with  Intel  Core  i7-2600  CPU  (8-core,  3.40GHz) 
and  another  with  AMD  Opteron  Processor  6134  (32-core,  2.30GHz).  Both  of  them  are  running 
Ubuntu  12.04  LTS  in  which  uses  EGLIBC-2.15  for  the  C  standard  library  implementation.  We  use 
MPFR-3.1.1  and  g++-4.8.1  C++  compiler  in  the  experiments. 

The  experimental  results  are  summarized  in  table  2. 

1 .  The  implementations  of  sin,  cos,  tan,  cosh,  sinh,  tanh  and  exp  functions  have  severe  prob¬ 
lems  when  used  with  non-default  rounding  modes  (toward  oo,  toward  —  oo,  and  toward  0). 
It  is  also  not  rare  to  have  the  problematic  case  in  practice.  For  instance,  cosh1  function  gives 
wrong  results  almost  50%  of  cases  (49,842  out  of  100,000). 

2.  We  have  not  observed  any  problem  under  the  default  rounding  mode  (toward  nearest  rep¬ 
resentable  number).  Also  the  implementations  of  acos,  asin,  atan,  tanh,  log,  logio,  and  sqrt 
functions  pass  our  tests. 

3  Patch  in  EGLIBC-2.17 


101 

102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 
113 


In  EGLIBC-2.17,  they  provided  a  patch  for  the  problem.  The  following  is  a  part  of  the  new 
implementation  of  sine  function  (IEEE754  double-precision)2: 

eglibc-2 . 17/libc/ sysdeps/ ieee754/dbl-64/s_sin  .  c  - 


_ sin (double  x){ 

double  xx,  res,  t,  cor,  y,  s,  c,  sn,  ssn,  cs,  ccs,  xn,  a,  da,  db,  eps ,  xnl ,  xn2 ; 

#if  0 

double  w [ 2 ] ; 

#endif 

mynumber  u,v; 
int4  k, m, n; 

#if  0 

int4  nn; 

#endif 

double  retval  =  0; 

S  E  T_RE  S  TORE_ROUND_5  3  B I T  ( FE_TONE ARE S T ) ; 


We  find  that  the  patch  does  not  really  fix  the  problem.  At  line  1 13,  it  resets  the  rounding  mode 
to  “round  to  nearest”  and  compute  the  value  of  sin(x)  while  ignoring  the  user-specified  rounding 
mode.  We  found  a  case  in  which  the  value  of  sin^,(x)  is  smaller  than  the  value  of  sinMPFR(x), 
which  violates  the  semantics  of  “toward  +oo”  rounding  mode: 

sinj,(— 3.93799)  =  0.714841448083829766879659928236 
sinMPFR  (-3.93799)  =  0.714841448083829771665831705916 

'Source  code  is  available  at  https://github.com/soonhokong/fp-test 

2 Available  at  http :  / / www .  eglibc  .  org/cgi-bin/viewvc .  cgi /branches/ eglibc-2_17/ 

libc/ sysdeps/ ieee754/ dbl- 64 /s_sin . c?view=markup 
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4  Conclusion 


We  report  serious  bugs  in  floating-point  computations  for  evaluating  elementary  functions  in  the 
Embedded  GNU  C  Library.  It  is  not  a  negligible  numerical  error  but  either  a  significant  error  (220 
ULP)  or  mathematically  incorrect  result  (i.e.  sin(x)  >  1053)  which  can  trigger  severe  problems  in 
the  following  computations.  Moreover,  the  chances  of  having  these  results  are  not  rare  at  all  as  we 
have  shown  in  section  2.  The  current  fix  does  not  mitigate  the  problem  but  hides  it. 
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